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A general expression for the Green's function of a system of N particles (bosons/fermions) 
interacting by contact potentials, including impurities with Dirac-delta type potentials is derived. 
, In one dimension for N > 2 bosons from our spectral determinant method the numerically calculated 

energy levels agree very well with those obtained from the exact Bethe ansatz solutions while they 
are an order of magnitude more accurate than those found by direct diagonalization. For N — 2 
bosons the agreement is shown analytically. In the case of N = 2 interacting bosons and one 
impurity, the energy levels are calculated numerically from the spectral determinant of the system. 
The spectral determinant method is applied to an interacting fermion system including an impurity 
to calculate the persistent current at the presence of magnetic field. 
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I. INTRODUCTION 



Due to the remarkable progress in the understanding of many features of mesoscopic physics the role of electron- 
electron interactions has been actively studied in such systems. Two frequently studied problems are the interaction 
induced derealization of two particles in disordered systems (we shall use the abbrevation TIP as it is common in the 
literature for two interacting particles) (I|E| and the persistent current (PC) first predicted by Buttiker et al. j3j and 
I | observed in metallic rings by Levy et al. and Chandrasekhar et al. fa]. Commonly used numerical approaches to 
these problems are the Green's function (GF) method ||U, the diagonalization method P-[TT||, the decimation method 
Jl2| and density-matrix renormalization-group (DMRG) method Jl3| . The many-body systems including impurities 
are frequently modeled by the Anderson-Hubbard tight-binding Hamiltonian [^) given by interacting electrons (with 
on-site interaction), and random potentials (with random site energy) of the impurities on a lattice. In the TIP 
problem the GF of the full Hamiltonian is obtained from Dyson's equation considering the interaction Hamiltonian as 
■ a perturbation. The two-particle localization length is then defined via a particular trace of the GF. In the PC problem 
the interplay between the interactions and impurities has generated much interest and several theoretical methods 
£^ are developed (see e.g. Ref. [|l4j and references therein) . To find the persistent current one needs to calculate the flux 
| dependence of the ground state energy which can be obtained from the lowest pole of the GF as a possible approach. 
t— ( i It is therefore desirable to construct the GF of the system which includes both impurities and the electron-electron 
interaction in order to treat the TIP and PC problems in a common framework. 

In this paper we develop a formalism to determine the GF of N-particle systems in which both the potential of 
the impurities and the interaction between the particles are given by contact potentials, i.e. by Dirac delta functions. 
Unlike in the case of the Anderson tight-binding Hamiltonian, in our calculation the positions of the electrons are not 
restricted to lattice sites but can be continuous coordinates. In our general formalism both the impurities and the 
interactions are treated as perturbations where the interaction or impurity strengths are the small parameters, and 
the GF is obtained from Dyson's equation. The resulting series for the GF can be expressed in terms of the GFs of 
the noninteracting many particle systems without impurities. If both the impurities and the interaction are given by 
contact potentials, the series can be given in a closed form, namely as the ratio of two determinants including only 
the GF of the unperturbed system. A similar (i.e. ratio of two determinants) form of the GF was found by Grosche 
, [15j for a system containing only one electron interacting with impurities. We have used this form of the GF in Ref. 
[16 1 to investigate the diffraction of the electron by the impurities, and in Ref. to study the crossover behavior 
of the localization problem in a waveguide with periodically placed identical point- like impurities. In this paper an 



extension of Grosche's approach to treat many-body interacting systems including impurities is presented. In our 
formalism we consider both the spinless bosons and the spin-full fermions. The many-body wave functions are the 
symmetrized products of the one-particle eigenstates of the system. In our calculations the bosonic/fermionic feature 
of the many-body system is fully incorporated. The energy eigenvalues of many-body systems are given by the poles 
of the GF. It implies that the energy levels are the zeros of the determinant in the denominator of the GF. We claim 
that the energy levels determined through this spectral determinant are much more accurate than those found by the 
traditional diagonalization methods. Therefore, here we propose the method of spectral determinants (SD) as a new 
approach to this problem. As it is seen below the indices of the matrix elements of the determinants are continuous 
variables. To evaluate the necessary SD in this case, one can use any complete orthogonal set as a basis. In our 
examples below this will be demonstrated using one-particle eigenstates. 
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As an application and to test our method, we consider the one-dimensional N-body problem in which the bosonic 
particles interact via a Dirac delta potential. It is well known that the exact result for the energy eigenvalues can be 
obtained from the Bethe ansatz ]T^-p2]|, so it is possible to compare our results with the exact ones. It is analytically 
proved that for the case of two interacting bosons the energy eigenvalues obtained from our method are identical 
to those found from the Bethe ansatz solutions. Moreover, in the case of three and four particles our numerical 
results agree very well with the Bethe ansatz solutions. Our results are compared with those obtained from direct 
diagonalization of the Hamiltonian, and it is demonstrated that in general our results are more than ten times more 
accurate. For higher energies it is shown that the direct diagonalization method is very inaccurate, which seriously 
limitates its applicability in numerical calculations. 

To see how much more effective our method is, we shall consider two nontrivial examples as a further application. 
We present our results for the system of two interacting bosons and one impurity, and for three interacting spin-full 
fcrmions and one impurity. In these cases no Bethe ansatz type of solution is known. Increasing the number of 
one-particle wave functions the energy eigenvalues converge. This way all the energy levels are calculated to four 
significant figures, and are taken as the exact energy levels. The errors of the energy eigenvalues found from our 
SD method and from the direct diagonalization of the Hamiltonian are compared using the same number of one- 
particle wave functions. In general, our method gives again an order of magnitude more accurate results than the 
diagonalization of the Hamiltonian. The accurate energy levels are necessary for calculating e.g. the persistent current 
which is the derivative of the ground state energy with respect to the applied magnetic flux enclosed by the system. As 
a demonstration we present our calculation of the persistent current for the interacting three-particle fermion system. 

Since in the above mentioned TIP and PC problems the Anderson-Hubbard tight-binding Hamiltonian can be 
regarded as a discretized version of the continuous model with Dirac delta potentials for both the interaction among 
the particles and with the impurities, we expect that our SD method might be applied successfully in numerical 
simulations. The work along this line is in progress. 

The rest of the paper is organized as follows. In section || we give ou r mo st general form of the GF for interacting 
N-particles systems including a finite number of impurities. In section III the GF is given for spin-full case. In all 
these sections our results for the GF are given as the ratio of two determinants expressed in terms of the GF of 
the noninteracting N-particle system. In section [Iv| our method is tested for interacting bosons for which the exact 
Bethe ansatz solutions are known. Our results are compared with the exact ones obtained from the Bethe ansatz 
solution and with those found from the numerical diagonalization of the corresponding Hamiltonian. In section [y] we 
calculate the energy levels for a system of two interacting bosons plus one impurity and compare the results with those 
obtained from the direct diagonalization of the Hamiltonian. In section VI the calculation of the per sistent current is 
presented for three interacting fermions and one impurity. Our conclusions are given in section VII , In Appendix ^ 
the derivation of the GF is given for the case of two interacting particles and M impurities. In Appendix |B we show 
that from our method the exact Bethe ansatz solution can be obtained analytically for two interacting bosons. 



II. MANY-BODY BOSON SYSTEMS WITH INTERACTION AND IMPURITIES 

In this section we consider the most general case, namely the system of ./V interacting spinless bosons plus M 
impurities with different strengths of potentials. To highlight the method to be applied it is useful to discuss the 
simple system of two noninteracting bosons with one impurity. It turns out that the full GF is a ratio of two struc- 
tured determinants. However, both determinants are so called 'continuous matrices' i.e. their indices are continuous 
variables. This feature of the matrices arises from the many body character of the system. 

The GF of the Hamiltonian H = Hq + Hi is given by 

G=(E-H -H 1 y\ (1) 
The Hamiltonian Hq of N identical noninteracting particles can be written as 

N 

ffo(x) = X>(*i); (2) 

i=l 

where the position of the zth particle is denoted by Xi, and x = (xi, x%, . . . , Xn). In general, the Xi are d dimensional 
vectors. In the examples presented below only one-dimensional systems are studied. The one-particle Hamiltonian 
h(x) is given by 

h(x) = -V 2 + V(x), (3) 
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where Vix) is the potential in the one-particle problem. Hereafter we shall use h = 2m = 1 units, where m is the 
mass of the particles. We assume that both the eigenstates and the energy eigenvalues are known for the Hamiltonian 
h of the one-particle problem. In many applications the potential V(x) is taken to be zero inside the box and infinity 
at its boundaries. The other typical case is when the potential V(x) is zero inside the box and periodic boundary 
conditions are applied. In this paper the latter is used in our examples. 

The Hamiltonian H± represents the effect of the impurity, and can be written as 

N 

#i(x) = k^2S(x p -u), (4) 
p=l 

where the position of the impurity is denoted by it, while k is the strength of the Dirac delta potential. The full GF 
of Hq + Hi satisfies Dyson's equation 

G = Go + G HiG, (5) 

where Go is the GF of the Hamiltonian Hq . 

To proceed further, it is instructive to consider only two noninteracting particles, i.e. N = 2. The result can be 
easily generalized then to the case N > 2. For two noninteracting particles the GF in coordinate representation is 
given |23| by 

n l \ I i \ ST ' i l } n{xi,X2)il)* rl {x'nx' 2 ) 
G {xi,x 2 \x 1 ,x 2 ) = E-E ' 

n 

where ip n (xi, x 2 ) is the nth eigenstate of the two noninteracting particles, and the corresponding eigenvalue is E n . In 
case of bosons the wave function is symmetric for the permutation of its variables, i.e. ip n (x±, x 2 ) = ip n (x 2 , x{). 
The iterative solution of Dyson's equation results in the series 



G(xi,x 2 \xi,x 2 ) = G (xi,x 2 \x 1 ,x 2 ) + / G (xi,x 2 \yi,y 2 ) Hi(yi,y 2 ) G (yx, V2\x 1 ,x 2 ) dyidy 2 

+ / G (xi,x 2 \yi,y 2 ) Hi(yi,y 2 ) G Q (yi,y 2 \y 3 ,y A ) Hi(y s ,yi) G {y3,yi\x'- l ,x' 2 ) dyidy 2 dy 3 dy 4 H . (7) 



After inserting Hi given in Eq. (|]) into Eq. (Q), one can perform some of the integrals in the series due to the 
presence of Dirac delta potentials. Finally, making use of the fact that for bosons 

G(xi,x 2 \x'i,x' 2 ) = G{x 2 ,xi\x' 2 ,x'i) = G(x 2 ,xi\x'i,x' 2 ) = G(xi,x 2 \x' 2 ,x[), (8) 

the series in Eq. (0) leads to 



G(xi,x 2 \x[,x 2 ) = G (xx,x 2 \x'i,x 2 ) + 2k J G Q (xi,x 2 \y,u) G Q (y, u\x\, x 2 ) dy 

+ (2k) 2 / G {xi,x 2 \yi,u) G (yi,u\y 2 ,u) G (y2,u\x'i,x 2 ) dyidy 2 -\ . (9) 



This is a Neumann series ]24j , known from the theory of integral equations. It can be rewritten in the more compact 
form 

G(xi,x 2 \x'i,x 2 ) = G (xi,x 2 \x'i,x 2 ) - J G (xi,x 2 \yi,u) K^ 1 (yi\y 2 ) G (y 2 ,u\x' 1 ,x 2 ) dy x dy 2 , (10) 
where the operator K(yi\y 2 ) is given by 

K{yi\y 2 ) = -77- ${yi - y 2 ) + G {yi,u\y 2 ,u). (11) 

For calculating the inverse of K one can choose a suitable complete orthogonal set to represent the operator if as a 
matrix. Using the identity for hypermatrices 



det ^ ° h d) = det( ^ det (° ~ bd ~ lc ) 



(12) 
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(which holds for arbitrary square matrices a and d, with detd 7^ 0), the full GF can be written as the ratio of two 
determinants 



dct 



G(x 1 ,x 2 \x' 1 ,x' 2 ) 



Go{xi,x 2 \x' 1 ,x' 2 ) 
G (y 1 ,u\x' 1 ,x 2 ) 



G (xi,x 2 \y 2 ,u) 
K(yi\y 2 ) 



detK(yi\y 2 ) 



(13) 



where the matrix in the numerator is discrete in its first row and column, but the rest is indexed by the continuous 
variables yi,y 2 . The 'continuous matrix' K can be given in matrix representation. Examples for treating this type of 
matrices will be given in Sees. |v|, [v| and [v|. Rewriting equation (10) m the form given by Eq. ( |l3|) is formal since 
in the numerator of Eq. (13) the dimension of the 1, 1 matrix element of the determinant is unity. However, later it 
turns out that in the more general case in which more than two particles are included as well as many impurities, the 
GF can be given similarly as a ratio of two determinants. 

Turning to the more general case of N > 2, the necessary steps to obtain the full GF are similar to those made in 
case of N = 2. The crucial point is to make use of the identity of the bosonic GF 



G(x|x') = G(Px|P'x'), 
where Px is an arbitrary permutation of the set (x\ 1 x 2l ...,Xn). One then finds 



(14) 



det 



G(x|x') = 



G (x|x') 
G («,y|x') 



Go(xKy') 
K(y\y') 



dettf(y|y') 



where 



K(y\y') = ~ - y') + G (u, y\u, /), 

Nk 



(15) 



(16) 



and y = (j/i, y 2) •••> VN-i)- The structure of the GF is similar to that given in Eq. (113]) e xcept that now the matrix 
K depends on some multidimensional variables. The indices of the 'matrix' K in Eq. (jig) are 'continuous' variables. 
The appearance of this kind of infinite-sized matrix is a consequence of the many-body feature of the Hamiltonian. 

We now turn to the case of interacting bosons. The previously presented method is suitable to handle the more 
complicated systems in which the interactions between particles are also taken into account. The derivation of the 
full GF of the two interacting bosons with many impurities is very similar to that shown before for the case of 
noninteracting particles, and the details are given in Appendix [a|. 

In the most general case the Hamiltonian Hq of the noninteracting N-particle system without impurities is given 
by Eq. (j|) . The rest of the total Hamiltonian of the system with M impurities and Dirac delta interactions with N 
particles is 



#i(x) 



M N N 

Ki > S(x p - ui) + A V S(x p - x q 



i=l p=l 



p<q 



(17) 



where x p are the positions of the bosons p = 1, • • • , N, Ui and Ki is the position and the strength of the impurity 
i = 1, • • • , M , and A is the strength of the interactions between particles. For simplifying the notations we shall use 
the vector x = (x\, x 2 , xn)- 

To find the full GF, a similar procedure given in Appendix |a] can also be carried out in this case, and the final 
result is 



det 



G(x|x') 



G (x|x') G (x|yi,yi,y') G (x| Ml ,y') 

Go(yi,yi,y|x') L{yi,y\y' x ,f) G {yi, yi, y>i, y') 

GoK,y|x') GoK,y|^,^,y') #i(y|/) 

G (u M , y|x') G Q (u M , y|yi,yi,y') G (wm, yK,y') 



Go(x|u M ,y') 
G {yi,yi,y\u M ,y') 
G (ui,y|u M ,y') 

K M (y\y') 



det 



L (yi,y\y'i,y') Go(?/i,?/i,yK, y ') 
Go(«i,y|yi,yi,y) #i(y|y') 

G (uM,y|yi,2/i,y') G (uM,yK,y') 



G (yi,?/i,y|wM,y') 
G (wi,y|«M,y') 

K M {y\y') 



(18) 
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where 



tfi(y|y') = -TT- ^^(y-yO + GoKyK,/), (19) 
L( yi , M, f) = - N(N 2 _ l)x t(yi v'i) t {N ~ 2) (f f) + Go(yi,yi,y\y[,y[,n (20) 

The variables y = (j/i , J/s, 2/4, yjv) an d y = (y3,Vi, 2/jv) (and the corresponding primed ones) appearing in the 
operators L and are (TV — l)-component and (N — 2)-component vectors, respectively. 

At first sight the final form of the Green's function seems to be very complicated but the matrices have a simple 
structure. The GF is again expressed as the ratio of two determinants. The variables of the full GF appear only in 
the first row and first column of the matrix in the numerator. Omitting the first row and first column of the matrix 
in the numerator the remaining matrix is the same as that in the denominator. The matrix in the denominator has 
an obvious structure, too: the 1,1 element of the matrix, L contains only the strength of the interaction between the 
two particles, while omitting the first row and the first column of this matrix it contains only the diagonal elements 
of the remaining matrix, i.e. the Ki contain the strength of the impurities. The other matrix elements depend only 
on the Green's function Go of two noninteracting particles. 

In summary of this section, the energy levels of the interacting boson systems with impurities are the poles of the 
Green function of the systems which coincide with the zeros of the denominator in Eq. ([l8|). This kind of form of 
the GF is a consequence of the contact potential assumed for the interactions between particles and the potential of 
the impurities. The denominator contains the operator K and L related to the impurities and the interactions of the 
bosons, respectively. In the special case when only the noninteracting N > 2 bosons with impurities are considered, the 
denominator reduces to det K (y|y') where K is given by Eq. (|l6|). The explicit form of these operators always available 
by using the symmetrized eigenfunctions of the noninteracting system. Then, these operators can be represented by 
an infinite dimension matrices. Note that, any basis set can be used (which is not necessaraly a symmetrized one) for 
evaluating the matrix elements of K and L. However, for numerical calculations of the matrices K and L one has to 
use a finite basis set. We would like to stress that in our method the truncation of the infinite basis set is the only 
approximation that is made. The advantage of our method against the direct diagonalization of the full Hamiltonian 
is that in our approach the bosonic/fermionic feature of the many-body systems is already taken into account through 
the Green function of the Hamiltonian Ho of the noninteracting, impurity free system. Then, it is ensured that the 
Green function of system of the interacting particles with impurities is also symmetrized according the nature of the 
particles. As we shall see below, the truncation introduce less error in the energy levels compared to that obtained 
by direct diagonalization than the lack of this symmetrization. 



III. SPIN-FULL FERMIONS WITH INTERACTION AND IMPURITIES 



In this section our spectral determinant method is generalized to an interacting many-body systems in which the 
particles can be bosons or fermions with spins. The Green function can be derived in a way similar to the spinless 
bosons discussed in the previous sections. Note that for spinless fermions with the Dirac-delta interaction potential 
the energy levels are the same as in the case of noninteracting fermions. 

The GF now depends on the spin indices and for noninteracting particles GF is 

where s = (si, S2, ■ ■■Sn) is a compact notation of the spins of the N particles. The interacting part of the Hamiltonian 
again consist of two terms: 

M N N 

#i(x) s | s / = ^2 K i y^J( x p - u i)$ s ,s> + A ^2 S(x p - x q )S StS t, (22) 

i— 1 p— 1 p>q=i 

p<q 

where all the terms are diagonal in the spin indices. Here ui and Ki are positions and the strength of the impurity 
i = 1, • ■ • , M, and A is the strength of the interactions between particles. 

Following the method developed in the previous sections we find that the GF for the system of interacting particles 
(fcrmion or bosons) with spins and impurities is given by 
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det 



G(x|x') sk 



G (x|x')s|s' G (x|?/i,y / 1 ,y') s |g, 

G {yi, y|x')s| S ' L(y%, ybi,y')sis' 



G (x|u 3 ,y') s |s' 
Go(2/i,2/i,yK-,y')s|i' 



. G (ui,y|x')s|s' GoCu^ylyi.y^yOsis' ^i,i(y|y / )s|s' 


det 


L(yi, y|yi, y')§is' G (yi, yi, y|%, y')s|s' 
GoK,y|yi,yi,y)s|s' *Q,j(y|y')s|s' 





where 



and 



SuS^-Vfr - y') C? + Go(«i,y|tii,y')i|r, 



>(yi,yVi,yV = " maT^iu ^ 1 -^ * (iV ~ 2) (y-y') C' + £o(yi,yi,y|y'i,yi,y')s| 



(23) 



(24) 



(25) 



Here the same notations are used for y = (2/1,2/3,2/4, ■■•,?/A r ) an d y = (2/3,2/4, ■••,2/iv) as m Eq. (|18|). The s,s denote 
the AT spin indices. 



IV. THE TEST OF THE SPECTRAL DETERMINANT METHOD FOR BOSONS 



In this section we apply our formalism in one dimension using periodic boundary conditions to two, three and 
four interacting spinless bosons without including impurities, as well as two bosons with one impurity. The energy 
eigenvalues of these many-body systems can be found from the zeros of a SD derived from GF formalism. In these 
systems one can calculate the exact energy levels from the well known Bethe ansatz solutions. Thus, our results 
can be compared with the exact ones. Using one-particle wave functions one can build up the symmetrized many-body 
wave functions, and in this truncated basis the Hamiltonian of the system can be diagonalized numerically. We use 
the same number of one-particle wave functions in the SD method as in the diagonalization of the Hamiltonian. If 
the number of one-particle wave functions used in the SD method is p, then for N > 2 interacting particles one needs 
to find the zeros of a SD of size p 1 . On the other hand, in the diagonalization of the Hamiltonian using the same 

p + N-l 



p one-particle wave functions the size of the matrix that should be diagonalized is 



N 



It is interesting to mention that for two interacting particles the energy levels are exactly the same as those obtained 
from the Bethe ansatz solutions. The details of the derivation is given in Appendix pi In sections IV A and IV B 



the derivation of the equations for the energy eigenvalues and the numerical results are presented for three and four 
interacting particles. The comparison of the SD method with direct diagonalization shows that the SD method leads 
to an order of magnitude improvement in the accuracy even in the worst case. 



A. Three interacting bosons in one dimension 



Now the method developed in section |n] is applied for A^ = 3 interacting particles in a one-dimensional box with 
periodic boundary conditions. To determine the Green's function Go for the case of three particles, it is again 
convenient to use the one-particle wave functions and the corresponding eigenvalues given by Eqs. (Bl) and 
Hereafter we shall use the length unit a = 1 for the size of the box. Then the Green's function Go is given by 



Go(xi,X2,x 3 \x' 1) x' 2 ,x' 3 ) = ^2 

{n 1 n 2 n 3 } 



VViin 2 n 3 {Xl, X2, X3) ^* ni n 2 n 3 ( x 'l , X 2i x i) 

E E ni E TL2 E n3 



(26) 



where the VVnn 2 n 3 ( x i, x 2, £3) are the eigenfunctions of the system of three noninteracting particles, and ni, 712,713 = 
0, ±1, ±2, • • •. For bosons these wavefunctions can be built up from the one-particle wave functions given by Eq. (Bl) 
Finally, the GF can be simplified as 



G (x 1 ,X2,x 3 \x' 1 ,x 2 ,x' 3 ) = — 



lpni{xi) 4>n 2 (x 2 ) i>n s {x 3 ) C P(1) (^'l) VC (2) (4) C P(3) (4) 



n-i_n 2 n 3 PeS 3 



En 3 



(27) 



where S3 denotes the permutation group of the numbers 1,2,3, and P is an element of the group, while P{i) is the 
number in the ith place in permutation P. 



G 



The energy levels of the system are the poles of the GF given by Eq. ([18]). From the form of the GF we can see 
that obtaining the poles is equivalent to finding the solutions of the equation det L — for E, where L is given by 
Eq. ©, i.e.: 

(28) 



det L(y, x\y', x') = det ( -— S(y - y') S(x - x') + G (y, y, x\y', y', x'; E) ) = 0, 



where A is the strength of the intercation between bosons. Using Eqs. (Bl) and ( |27| ) the second term in det L becomes 

1 f ,2-KiUn 1 +n2)(.y-y')+n 3 (x-x')] , n 2Tri\n 1 (y-x')+n 2 (x-y')+n 3 (y-y')] 

G (y,y,x\y',y',x') = ^ V - tt> B ■ (29) 



E E ni E n2 E n3 



The determinant of the operator L in Eq. (|28|) can be rewritten as a determinant of an infinite sized matrix using 
any convenient basis (which is not necessarily that of symmetrized wave functions). In our calculation wc chose 



/ \ 2-Kiiky+lx) 

(fkl{y,x) — e y ST ; . 
Thus, the operator L in this basis can be given by 

L(kl\k'l') = J <p* kl (y,x) L(y,x\y',x') ipk'l'(y', x') dy dx dy' dx' . 

After some tedious but straightforward algebra one can find that 



1 



1 



L(kl\k'l') = -— S kX S ltV + - 



2 Sk+i. 



k'+V 



E-E,, 



Ei — Ek-i' 



E 



1.1' 



E — Eh-n — E n — El 



(30) 



(31) 



(32) 



The sum in this equation has a similar structure to Eq. (BIO), therefore using the identities given in Eq. (Bll) the 
summation can be carried out yielding 



L(k,l\k',l') 



S k ,k'Sl,l> (-ax + T^FI cot I z ) + A ( E '> fc ' for k even ' 

S k ,k'Si,i' (31 + tan f z) + A(E; k, l\k', I'), for k odd, 



(33) 



where 



A(E;k,l\k',l') = - 



?k+i,k'+i' 



3 E-4ir 2 [I 2 + I' 2 + (k - l') 2 Y 



E 
2^2 



2l 2 - k 2 . 



(34) 



One can see that choosing an appropriate basis it is possible to express all the matrix elements of the operator L 
in a closed form. We now solve Eq. (|28J) numerically to find the energy levels from our method and compare them 
with the Bethe ansatz solutions determined in Eqs. (B14) and ( B15| ). 

In numerical calculations we have to truncate the matrix L(k, l\k' , /'). To form the basis (fiki(x, y) given in ( [30] ) we 
chose 2m+ 1 one-particle states ipk(x) (k = 0, ±1, ±2, • • • , ±m) as given in Eq. (Bl). Therefore, we have (2m + 1) dif- 
ferent basis functions fki(x, y) with k, I = 0, ±1, ±2, • • • , ±m. Thus, the infinite sized matrix L(k, l\k', I 1 ) is truncated 
to a (2m + l) 2 x (2m + l) 2 matrix. 

The logarithm of the absolute value of det L(k, l\k' , I') is plotted as a function of the energy E in Fig. [j]. In this 
plot 2m +1 = 7 one-particle wave functions were used. The zeros of det L(k, l\k', I') (i.e., in the plot at the energy 
values where the logarithm of det L(k, l\k' , I') goes to —00) give the energy levels of the system. One can see from 
the figure that these energy values agree very well with those found from the Bethe ansatz solution. Note that the 
det L(k, l\k', I') has poles at the energy levels of the noninteracting system. This can be seen clearly from the figure. 
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FIG. 1. The logarithm of the absolute value of det L(k, l\k' , I') as a function of E for N = 3 interacting particles (solid 
line). The strength of the interactions is A = 4.0. The number of one-particle wave functions was 2m + 1 = 7. Including 
the interactions, the energy levels, i.e. the position of the zeros obtained from the Bethe ansatz solution are shown by dashed 
vertical lines. For clarity, we indicate the energy levels of noninteracting particles with dotted vertical lines (Ek — 4-K 2 k 2 , where 
fc = 0,l,2,.--). 

In Fig. H the relative errors of the first six energy levels obtained from Bethe ansatz solutions and our SD method 
using 2m + 1 one-particle wave functions are plotted as a function of m. One can see that the relative errors decrease 
with increasing m. It implies that including only a few number of one-particle wave functions in our SD method is 
enough to get a very satisfactory agreements with the exact energy levels obtained by Bethe ansatz solutions. 



1 

0.1 

g 0.01 
o 

0.001 
0.0001 
1e-05 



10.3284947 
57.1227798 
96.6011996 
101.131383 
128.763754 
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FIG. 2. The relative errors (in percentage) of the energy levels found (the different energy levels are indicated by different 
symbols) from the SD method using 2m + 1 one-particle wave functions as a function of m for the first six energy levels. The 
strength of the interactions is A = 4.0. 



In Fig. H we plotted the relative errors of the energy levels found from the SD method for the first 100 energy levels. 
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FIG. 3. The relative errors (in percentage) of the the first 100 energy levels found from the SD method using 2m + 1 = 19 
one-particle wave functions. The strength of the interactions is A = 4.0. 

Diagonalization of the Hamiltonian can be performed by using the same one-particle wave functions as in the SD 
method. In Table | we have listed the first few energy levels obtained from Bethe ansatz solutions, diagonalization 
of the Hamiltonian of the system and the SD method. It can be seen from Table || that the relative errors of the 
energy levels obtained from the SD method are an order of magnitude smaller than those of the diagonalization of 
the Hamiltonian. 

One can see that the errors of the first 50 energy levels are less than 0.01%, in spite of the fact that only 9 one- 
particle wave functions were used. It is surprising that only a relatively small number of one-particle wave functions 
needs to be used for obtaining a rather accurate result. The accuracy breaks down after the first 50 energy levels. 
This is the point where the unperturbed energy levels corresponding to the calculated levels start missing in case of 
9 one-particle base functions. 



Bethe ansatz Diagonalization Spectral determinant method 


E 


m = 2 


m — 3 


m = 4 


m = 2 


m = 3 


m = 4 


10.3284947 


10.7066 (3.7) 


10.5994 (2.6) 


10.5403 (2.05) 


10.3264 (0.02) 


10.3277 (0.008) 


10.3281 (0.004) 


57.1227798 


57.9505 (1.4) 


57.6827 (1.0) 


57.5488 (0.75) 


57.1508 (0.05) 


57.1196 (0.006) 


57.1508 (0.049) 


96.6011996 


97.4767 (0.9) 


97.1866 (0.6) 


97.0426 (0.46) 


96.5507 (0.05) 


96.6507 (0.051) 


96.6507 (0.051) 


101.1313830 


102.8964 (1.7) 


102.7003 (1.6) 


102.5885 (1.44) 


101.1061 (0.03) 


101.1231 (0.008) 


101.1277 (0.004) 


128.7637540 


129.3888 (0.5) 


129.1411 (0.3) 


129.0332 (0.21) 


128.7006 (0.05) 


128.7584 (0.004) 


128.7756 (0.009) 


137.2597820 


137.7723 (0.4) 


137.4139 (0.1) 


137.2811 (0.02) 


137.2505 (0.01) 


137.3005 (0.030) 


137.2630 (0.002) 


176.7382020 


177.5799 (0.5) 


177.0631 (0.2) 


176.8474 (0.06) 


176.6504 (0.05) 


176.7504 (0.007) 


176.7340 (0.002) 


219.5666420 


222.0227 (1.1) 


221.2907 (0.8) 


221.1095 (0.70) 


212.7003 (3.13) 


219.5003 (0.030) 


219.6003 (0.015) 


220.2550290 


222.4741 (1.0) 


221.8307 (0.7) 


221.6792 (0.65) 


220.1503 (0.05) 


220.2378 (0.008) 


220.3003 (0.021) 


254.5148790 


256.0261 (0.6) 


255.0061 (0.2) 


254.8096 (0.12) 


252.7502 (0.69) 


254.4002 (0.045) 


254.5502 (0.014) 



TABLE I. The energy levels from Bethe ansatz, diagonalization of the Hamiltonian and the SD method. In the diagonal- 
ization method we used the same one-particle wave functions as in the SD method. The percentage errors are indicated in 
brackets. The strength of the interactions is A = 4.0. 
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B. Four interacting bosons in one dimension 

In this subsection the equation for the energy levels is derived for TV = 4 interacting particles in a one-dimensional 
box with periodic boundary conditions. The necessary steps to find this equation are similar to those in the previous 
subsection. In fact, the SD method can be generalized straightforwardly to any number of interacting particles, 
although the algebra becomes more complicated for N > 4. The energy levels for N — 4 particles are the zeros of the 
determinant of operator L given in Eq. ( 

The first few energy levels obtained from solving are listed in Table y. The strength of the interactions is A — 4.0. 
For the sake of comparison the exact energy levels obtained from the Bethe ansatz solution and the results found 
from the direct diagonalization of the Hamiltonian are also listed in the table. Both in the SD method and in the 
diagonalization method the same 2m + 1 one-particle wave functions were used. The percentage errors with respect 
to the exact Bethe ansatz solutions are indicated in brackets. 

From the table one can see that the relative errors using the SD method are more than an order of magnitude 
smaller than those obtained from direct diagonalization. In particular, using a rather small number of one-particle 
wave functions, 2m + 1 = 9, the ground state energy is much more accurate (0.01%) compared to those obtained from 
the diagonalization method. 

Note that to find the energy levels with the above indicated errors only the determinant of a 9 x 9 matrix had to be 
calculated. This suggests that in numerical calculations the SD method is much more efficient than the conventional 
diagonalization method. Comparing the two methods one realizes that although in our method more algebra is 
needed than in the diagonalization methods, less numerical effort is required to obtain the values of the energy 
levels accurately. It is important to note that the expression for determining the energy eigenvalues becomes more 
complicated with increasing the number of interacting particles but the necessary steps are quite straightforward and 
a quite simple algorithm can be given. 



20). 



Bethe ansatz 



Diagonalization 



Spectral determinant method 



E 


m = 2 


m — 3 


m = 4 


m = 2 


m — 3 


m — 4 


20.8016 


0.00 (100.0) 


21.34 (2.6) 


21.22 (2.01) 


20.7904 (0.05) 


20.7976 (0.02) 


20.7997 (0.009) 


70.8207 


72.18 (1.9) 


71.74 (1.3) 


71.51 (0.97) 


70.7561 (0.09) 


70.8031 (0.02) 


70.8135 (0.010) 


113.7688 


115.36 (1.4) 


114.83 (0.9) 


114.56 (0.69) 


113.6250 (0.13) 


113.7380 (0.03) 


113.7576 (0.010) 


118.4289 


120.12 (1.4) 


119.59 (1.0) 


119.31 (0.74) 


118.3495 (0.07) 


118.4014 (0.02) 


118.4168 (0.010) 


149.7776 


151.31 (1.0) 


150.78 (0.7) 


150.51 (0.49) 


149.4374 (0.23) 


149.7434 (0.02) 


149.7660 (0.008) 


158.5862 


160.32 (1.1) 


159.78 (0.8) 


159.49 (0.57) 


158.3920 (0.12) 


158.5536 (0.02) 


158.5722 (0.009) 


178.7153 


179.98 (0.7) 


179.47 (0.4) 


179.25 (0.30) 


178.3744 (0.19) 


178.6877 (0.02) 


178.7084 (0.004) 


191.2125 


193.25 (1.1) 


192.43 (0.6) 


192.08 (0.46) 


191.0778 (0.07) 


191.1803 (0.02) 


191.2018 (0.006) 


195.3592 


197.32 (1.0) 


196.46 (0.6) 


196.20 (0.43) 


195.2056 (0.08) 


195.3277 (0.02) 


195.3468 (0.006) 


237.5430 


240.82 (1.4) 


239.00 (0.6) 


238.59 (0.44) 


237.7028 (0.07) 


237.3591 (0.08) 


237.5022 (0.017) 


238.2730 


241.38 (1.3) 


239.68 (0.6) 


239.31 (0.43) 




238.2151 (0.02) 


238.2512 (0.009) 


276.3426 


279.69 (1.2) 


277.84 (0.5) 


277.40 (0.38) 


274.9712 (0.50) 


275.8985 (0.16) 


276.2966 (0.017) 


277.9374 


280.87 (1.1) 


279.30 (0.5) 


278.97 (0.37) 


277.7805 (0.06) 


277.8731 (0.02) 


277.9146 (0.008) 


281.8038 


281.29 (0.2) 


283.42 (0.6) 


282.94 (0.40) 


285.7567 (1.40) 


281.5560 (0.09) 


281.7657 (0.014) 



TABLE II. The exact energy levels obtained from the Bethe ansatz and our numerical results using 2m + 1 one-particle wave 
functions. In the diagonalization method we used the same one-particle wave functions as in the SD method. The relative 
errors in percentage are indicated in brackets. The strength of the interactions is A = 4.0. 
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V. TWO INTERACTING BOSONS WITH ONE IMPURITY IN ONE DIMENSION 



As a further illustration of the SD method, in this subsection we consider a system including two particles and an 
impurity in a one-dimensional box with periodic boundary conditions. In this case no exact results, such as Bethe 
ansatz solutions, are known for the energy levels. The Hamiltonian of the system is H = Hq + Hi, where Hq is given 
by Eq. (0). The interaction Hamiltonian of the system can be written as 



H 1 (xi,x 2 ) — XS(xi — x 2 ) + k[8 {x\ — u) + 8 (x 2 — u)] 



(35) 



where the strength of the interaction between the two particles and the strength of the potential for the impurity are 
denoted by A and k, respectively. The impurity is located at u. The GF of the system is given by Eq. (fL8|) with 
M = 1. The energy levels of the system are the roots of the denominator in (jl^), which in our case, can be written as 



det 



L(y,y';E) G (y,y\y> ,u; E) 
G (y,u\y',y';E) K(y,y';E) 



0. 



where 



(36) 



L(y,y') 



■t8(v-v') +G (y,y\y',y'), 



K(y,y') = ~TT S (y ~ v') + G (y,u\y\u). 
Zk 

The Green's function Gq{x\, x 2 \xi, x' 2 ) of the Hamiltonian Hq can be written as 

VVpEl) ?Pn 2 (x 2 ) l}}* np{i) {x' 1 ) 1p* npm {x' 2 ) 



Go(xux 2 \x[,x' 2 ) = y 

n 1 n 2 PeS 2 

= 1 (y- tpnjxi) 1pm{x 2 ) lpn( x l) V4(4) 

2 If- E — E n — E m 



E — E ni — E n , 2 

tpn(xi) 1pm(x 2 ) Vm( x 'l) C04) 



E 



E — E n — E„ 



(37) 
(38) 



(39) 



where ipn{x) and E n are given by Eqs. (Bl) and (B2), respectively (with a — 1). In Eq. (36) the four matrix elements 
contain the two-particle Green's function Gq(xi, x 2 \x' x , x' 2 ) with different arguments. The four matrix elements are 
operators but for numerical calculations we need to represent them by matrices. This can be done by using a set of 
orthogonal wave functions: 



Vk{y) 



— cZmkiy-u) 



(40) 



where the extra phase e 2Trm is introduced for the sake of convenience. Then, in this basis, the Green's function Go 
appearing in the matrix elements 1, 1, 1, 2, 2, 1 and 2, 2 of Eq. ( |36| ) becomes 



4 n) (fc|0 = 4,*]T 

GP(k\l) = 
G^(k\l) = 



1 



E — E n — Ek-n 



1 



GP(k,l)= 1 - 



E-Ei- E k -i 
1 

E — Ek — Ei-k 
1 



= G, 



(1,2) 



(l\k), 



(41) 



E — Ek — Ei 



Sk ' l ^2 E-E k -E^) 



The sum over n in G^ 1 ^ {k\l) can be carried out by using the identities Eq. (Bll), while the summation in G^ 1 (k\l) 
can be performed using the identity |25| 



(22), 



E 



7T cot 1TZ 



(42) 



Finally, the matrix representation of Eq. (pq) becomes 
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det 



where the four elements of the hypermatrix are 



Af(n) M < 12 ) 

M (21) M (22) 



M, 



M, 



(ii) 

k,l 

(12) _ 

E - 4tt 2 [P + (k - I) 2 } ' 
1 1 



1 



M, 



(21) 
fej 



M, 



(12) 



with 



2k 47TZ. 



■ COt 7TZ2 + 



S-47r 2 (fc 2 + P)' 



3 ^-cot(f2i) for fc even, 
~5rW tan (l z i) for k odd ' 



21 



J5 
2^2 



k 2 and Z2 



J5 
4^2 



k 2 . 



(43) 



(44) 



(45) 



(46) 



Note that to calculate numerically the determinant in Eq. (|43J) it is useful to apply an equivalent form of the general 
matrix identity for hypermatrices given in Eq. (n2|) 



det 



a b 
c d 



= det (a) det (d — ca 1 fy 



(47) 



(which holds for arbitrary square matrices a and d, with det a ^ 0). Then, the dimensions of the minors are half of 
the original one. Since the matrix L 1 ^^ is diagonal, it is easy to calculate its inverse. 



The first few energy levels obtained from solving Eq. (J43j) numerically are listed in Table |III| . The strength of the 
interactions and the strength of potential for the impurity arc A = 3.0 and n — 2.0, respectively. 



'Exact' E (4 digits) 


Diagonalization 


Spectral determinant method 


E 


m = 2 


m — 3 


m = 4 


m = 2 


m = 3 


m = 4 


6.1175 


6.3153 (3.2) 


6.2583 (2.3) 


6.2262 (1.78) 


6.1164 (0.018) 


6.1171 (0.0067) 


6.1173 (0.0032) 


46.8629 


47.1405 (0.6) 


47.0458 (0.4) 


47.0002 (0.29) 


46.8624 (0.001) 


46.8628 (0.0002) 


46.8629 (0.0000) 


50.4776 


50.8832 (0.8) 


50.7582 (0.6) 


50.6908 (0.42) 


50.4695 (0.016) 


50.4750 (0.0051) 


50.4765 (0.0023) 


82.9012 


83.0716 (0.2) 


82.9998 (0.1) 


82.9746 (0.09) 


82.9007 (0.001) 


82.9010 (0.0001) 


82.9011 (0.0001) 


85.5893 


85.8157 (0.3) 


85.7852 (0.2) 


85.7338 (0.17) 


85.5869 (0.003) 


85.5891 (0.0003) 


85.5893 (0.0001) 


91.2389 


91.7645 (0.6) 


91.5606 (0.4) 


91.4876 (0.27) 


91.2224 (0.018) 


91.2350 (0.0043) 


91.2373 (0.0018) 


165.5271 


165.9033 (0.2) 


165.7284 (0.1) 


165.6938 (0.10) 


165.5104 (0.010) 


165.5263 (0.0005) 


165.5269 (0.0001) 


169.3990 


170.0360 (0.4) 


169.7856 (0.2) 


169.6524 (0.15) 


169.3588 (0.024) 


169.3938 (0.0030) 


169.3972 (0.0011) 


203.1501 


203.5057 (0.2) 


203.3393 (0.1) 


203.2789 (0.06) 


201.3532 (0.885) 


203.1501 (0.0000) 


203.1501 (0.0000) 


206.9552 


207.5140 (0.3) 


207.2107 (0.1) 


207.1694 (0.10) 


207.1683 (0.103) 


206.9520 (0.0015) 


206.9547 (0.0002) 


207.1812 


207.5666 (0.2) 


207.4673 (0.1) 


207.3835 (0.10) 


209.2883 (1.017) 


207.1799 (0.0006) 


207.1811 (0.0001) 


210.9868 


211.9115 (0.4) 


211.4817 (0.2) 


211.3105 (0.15) 




210.9672 (0.0093) 


210.9816 (0.0024) 



TABLE III. The 'exact' energy levels (results to four significant figures using 61 one-particle wave functions) and our 
numerical results obtained from Eq. (j^j using 2m + 1 one-particle wave functions. In the diagonalization method we used the 
same one-particle wave functions as in the SD method. The percentage errors are indicated in brackets. The strength of the 
interactions and the strength of potential for the impurity are A = 3.0 and n — 2.0, respectively. 



No Bethe ansatz type solution is known when the system includes impurity. The 'exact' energy eigenvalues are 
calculated by increasing the number of one-particle wave functions until all the energy levels listed in the first column 
of the table are converged to four significant figures. To reach this convergence, 61 one-particle wave functions were 
used. Then the errors of the energy eigenvalues found from the SD method and from the direct diagonalization of the 
Hamiltonian are compared using the same number of one-particle wave functions. As it is seen from Table 
method usually gives again an order of magnitude more accurate results than the diagonalization of the Hami 



[II, our 



tonian. 
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VI. THREE INTERACTING FERMIONS WITH ONE IMPURITY IN ONE DIMENSION AND 
CALCULATION OF THE PERSISTENT CURRENT 



In this section our SD method is applied for a more complicated problem. We consider the one dimensional ring 
including N — 3 interacting 1/2 spin-fermions. The external magnetic field is applied along the axis perpendicular to 
the plane of the ring. In addition, a single impurity is also included in an arbitrary position. The persistent current 
is the derivative of the ground state with respect to the flux of the magnetic field inside the ring ||^6| : 

where E($>) is the energy of the ground state of the system with flux $ throu gh t he ring. The ground state energy is 



the smallest poles of the GF for fermion systems given by Eq. (|23| ) in section [III 
The GF is constructed from the one-particle wavefunctions: 

ipn,*(x)s = e 2 ™ (n+$) V s , (49) 

where a = |f) s = 1,2 denoting the spinor index of the spin state a and $ is the magnetic flux in units of 
flux quantum h/e. The corresponding eigenvalues of the free one-electron Hamiltonian are E n = 47r 2 (rt + <I>) 2 with 
n = 0,±l,±2,---. 

The unperturbed GF appearing in Eq. ( |23| ) can be constructed from the above given one-particle wave functions. 
For fermions the symmetrized form is: 



Il2,ff2 \ X 2) S^n 3 ,a 3 

G ° - 3! h. Jk. { } E-E ni - E n2 - E n3 



Pes 3 "i"2™3 



t7 l cr 2 cr 3 

(50) 



where S3 denotes the permutation group as given after Eq. ( p7| ) . The denominator of the GF in (|23|) can be transformed 
into a matrix form by using the non-symmetrized product of the one-particle wave functions as it was done in the 
case of bosons. The determinant of this matrix takes the following 2 by 2 hypermatrix form: 



det 



M (21) M (22) 



(51) 



where the p = (k, I, s\, S2, S3), q = (k' , I', s^, s' 2 , s' 3 ) matrix element of the matrices M is given by 



M (ll) _ _ J_r 1 r r r , 1 ^n 1 +n 2 ,k^n 3 ,lSn r(1) +np (2h k'Sn P(3h l' . , 

VI ~ 3^°fe,fe' '^' «i,«i°S2,« 2 s 3,4 + g Z^ Z^ E — E —E —E ° s P(.i)' s i s P(.2)' s 2 s P(3)t s 3^ \p£) 



P&S3 ni,n2,n3 



M (12) _ 1 S ni+n2! kSn 3 ,lS np(2h k'S np(3) M 

1V1 P1 ~ R Z^ Z^ W V) _ R _ E °SP(1),4^P(2),S 2 °SP(3):<,' \ 06 ) 

M (21) _ 1 \^ ^n 2 ,fc^n 3 ,i^n P(1) +np (2) ,fc' ^ P(3 ) ,i' - - „ , , 

P9 r Z^ Z^ E - E - E - E c, sp (1) ,s' 1 Osp (2) , S2 o S p ( 3 ) ,4, (.a*; 

PeS 3 ni,ri2,n 3 111 ™ 2 ™ 3 

,, r (22) 1 r X X X 1 1 ^n 2 ,fc<5n3,i^np (2) ,fc"5np (3) ,i' 

^ 3T d *.* ,d, ' , '0-i,»lO« a ,4 d «8.«i + « 2^ Z^ ~^~E~~E~~E V (1 ),<V( 3 ),« 2 ^i»(3).« 3 - i 55 ) 

PeS3m.na.n3 111 ™ 2 ™ 3 

Since these matrices are infinite dimensional ones we truncate them by restricting k, I between —m and to, where 2m+l 
is the number of one-particle wave functions given by (p)|). Because of the spin indices each (k,l), (k',l r ) elements 
of the matrices M still correspond to 8 by 8 matrices. Therefore in Eq. (|5l]) the hypermatrix is a 2 • 8 • (2m + l) 2 
dimensional matrix. Grouping these elements by the total S z of the three fermions we may reduce the dimension of 
the matrix into a (2 + 3) ■ (2m + l) 2 dimensional one. This dimensional reduction is possible because we omit the 
energy independent factors in the matrix elements. 

The above given matrix elements still contain summations which can be evaluated by using the following identity 
0: 
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OO 

E (n + j>)2 _ Z 2 = ^ ( cot *@ -*)- cot *@ + *))■ ( 56 ) 

71= — OO 

In Fig. H the persistent currents obtained by our SD method as a function of flux in units of flux quantum h/e are 
shown for different impurity strengths k. In this calculation 2m +1 = 5 one-particle wave functions are used. 




-0.5 -0.3 -0.1 0.1 0.3 0.5 

FIG. 4. The persistent current calculated by the SD method is plotted (in arbitrary units) as a function of the flux (in units 
of 4> = h/e) for fixed interaction strength A = 2 and different impurity strength k. We used 2m +1 = 5 one-particle wave 
functions. 



One can see that the persistent current suppresses with increasing the strength of the impurity. Our purpose in 
this paper is to demonstrate how our SD method can be applied to different problems such as the calculation of the 
persistent current. A more complete analysis of the persistent current problem is beyond the scope of this paper. The 
work along this line is in progress. 



VII. CONCLUSIONS 



In this paper the Green's functions of many-body systems were calculated when the interactions between the 
particles are Dirac-delta type contact potentials, and when the system contains impurities described by contact 
potentials. The Hamiltonian of the system is taken as a sum of two parts, the first of which contains only the sum 
of the one-particle Hamiltonians, while the second contains the rest, namely the interactions between the particles 
and the potentials of the impurities. Using the Dyson equation for the Green's function the iterative solution for 
the Green's function can be formally given as the ratio of two determinants. However, in this formal expression the 
matrix elements of the determinants are 'indexed' by continuous variables. To handle these determinants one can use 
any complete orthogonal set of wave functions and the corresponding operators in the determinants can be given in 
the usual matrix representation. This makes the numerical calculation of the Green's function possible. 

After presenting the general formalism for deriving the Green's function, several examples were discussed in order 
to show the relevant steps in the algebra as well as to test the method. Only the energy eigenvalues were calculated 
for different systems. However, in our method the full Green's function is given, therefore, in principle, it can be used 
for calculating other types of physical quantities. The energy levels can be found by determining the poles of the 
Green's function. Since the Green's function of the system is given by the ratio of two determinants, one only needs 
to calculate the zeros of the spectral determinant (in the denominator of the Green's function) to obtain the energy 
eigenvalues of the system. 

For testing our method, we used a one-dimensional model of bosons interacting via Dirac-delta interactions. In 
this case the energy eigenvalues can be calculated exactly from the well known Bethe ansatz solution . In the case 
of two interacting bosons it was shown algebraically that the SD method yields exactly the same energy levels as 
those obtained from the Bethe ansatz solution. For the case of three and four interacting bosons an equation for 
determining the energy levels was derived and then it was solved numerically. The results obtained from our method 
are in very good agreement with those obtained from the Bethe ansatz solution. A comparison of our method with 
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the direct diagonalization of the Hamiltonian was also presented. The same number of one-particle wave functions 
was used in the latter method as in the SD approach. It was demonstrated that our method gives, in general, ten 
times smaller relative errors than the diagonalization method. The energy of the ground state of the system can be 
obtained with even better accuracy in spite of the fact that only a few one-particle wave functions were used (typically 
9; the relative errors were then less than 0.01%). As a nontrivial example we calculated the energy levels of a system 
including two interacting bosons and one impurity. Both the interactions and the potential of the impurity were of 
Dirac-delta type. Applying our general method, we derived an equation for the energy levels of the system. Since 
no Bethe ansatz type of solution is known in this case, by increasing the number of one-particle wave functions we 
calculated the energies by making the result convergent to 4 significant figures. Taking these energy levels as the exact 
ones, we compared our method with direct diagonalization. A much better accuracy was found in the SD method 
(e.g. for the ground state with 9 one-particle wave functions the relative error was about 0.003%, while it was 2.6% 
with the diagonalization method). 

Our method was also applied to spin-full fermion systems. A better accuracy was achieved in calculating the ground 
state energy comparing the diagonalization method. Therefore our method is expected to be suited for numerical 
calculations for the persistent current. We demonstrated the applicability of our SD method for the case of three 
interacting fermions and one impurity, and calculated the persistent current. We defer a thorough study of the 
persistent current in a later publication. 

In summary, we developed a spectral determinant method with which much more accurate values of the energy 
levels can be obtained numerically than with the usual direct diagonalization method. In our formalism the full 
Green's function is available, thus it might be used in the TIP and persistent current problem mentioned in the 
introduction. 
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APPENDIX A: TWO INTERACTING PARTICLES WITH M IMPURITIES 



In this section an expression for the GF is derived for the case of two interacting particles and M impurities. The 
interaction Hamiltonian Hi is given by 

M 

Hi = Ki(s(xi - Ui) + S(x 2 - Uifj + A 5(xi - x 2 ), (Al) 

i=l 

where U{ is the location of the ith impurity, of strength K{ and A is the strength of the interaction between the 
particles. As in the case of a single impurity, the first-order correction (in Hi) of the iterative solution of the Dyson 
equation is 

= J G (xi,x 2 \yi,y 2 ) Hi(yi,y 2 ) G (yi,y 2 \x'i,x' 2 ) dyidy 2 . (A2) 

Inserting Hi into the above expression the integrals can be performed exactly for Dirac delta potentials. Collecting 
the identical terms arising from permutational symmetry, the first order correction of the GF is 

„ M 

= / ^2k 4 Go{xi,x 2 \y,Ui)G {y,Ui\x'iX 2 ) + A G (xi,x 2 \y,y)G (y,y\x'i,x 2 ) dy. (A3) 
J »=i 

Similarly, the second order correction takes the form 

G {2) = J G {xi,x 2 \yi,y 2 ) H 1 (y 1 ,y 2 ) G {yi,y 2 \y'i,y' 2 ) H 1 (y[,y / 2 ) G (y'i, y' 2 \x\, x' 2 ) dyidy 2 dy[dy' 2 . (A4) 
Higher order corrections can be easily found. In a way similar to the case of G^ we have 
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G (2) = j ^2 2k * G o{xi,x 2 \y 7 u i )2K j G (y,u i \y',Uj)G (y',Uj\a/ 1 X2) dydy' 

. M 

G (xi,x 2 \y,Ui)\ G (y,u i \y',y')G (y',y'\x' 1 x' 2 ) dydy' 



<=i 
M 



E A G , (a;i)»2|j/,J/)2Kj Go(y,2/|y',u J )Go(?/',u :) |x' 1 x2) dydy' 
A G (xi,x 2 |y,y)A G (y, y|y', y')G (y', y'lx'^) dydy'. 



(A5) 



It can be shown that the GF including the higher order corrections is 

, M+l „ M+l 



Mv) B a (y) dy+ Mv) C a p(y,y') B p {y') dydy' + . . 

a=l a,0=l 
. oo M+l 

= G °+ / E E A ^ { c <*?M)) n Mv') dydy', 



n=0 a, /3=1 



where 



A Q (y) = ( 2 ^G (x u x 2 \y,^)\ , = / G,(y' u 3 y r x' 2 )\ 

ayiJ! \ XG (x u x 2 \y,y) J 1 ^ ; ^ G (y', y'K> x' 2 ) J' 



are (M + l)-component vectors (i = 1, . . . , M) and 



r- f„ - ( 2k i G (y,u t \y' ,u 3 ) A G (y, Ui\y', y') 
c a p{y,y)- I 2k . GoMj/> .) a Go (y,y|y',y') 



is an (M + 1) x (M + 1) matrix. This is a Neumann series which can be summed: 

M+i 



G = Go+J A a (y) (6(y-y')6 a p-C af3 (y,y')) * B p (y') dydy'. 



(A6) 



(A7) 



(A8) 



(A9) 



a, 0=1 



Using Eqs. ( A7), ( A8), and the matrix identity in Eq. (|l2|), the GF can be rewritten as the ratio of two determinants: 



det 



G(x 1 ,x 2 \x' 1 ,x 2 ] 



Gq{xi,X2\x' 1 ,x' 2 ) G Q (xi,x 2 \y',y') G (xi,x 2 \ui,y') 

Go(y,yWi,^2) L (yW) G (y,y\ui,y') 

G Q (ui,y\x' x ,x' 2 ) G (ui,y\y',y') Ki(y\y') 

Go(uM,y\x'i,x 2 ) G (u M ,y\y',y') G Q (u M ,y\ui,y') 



G (xi,x 2 \uu,y') 
Go(y,y\uM,y') 

G (ui,y\u M ,y') 

K M (y\y') 



det 



L(y\y') G {y,y\ui,y') 
G (u u y\y',y') ^i(y|y') 

G {u M ,y\y',y') G {u M ,y\ui,y') 



Go(y,y\u M ,y') 
G (ui,y\uM,y') 

KmW) 



(A10) 



where 



and 



Ki{y\y') = -— 5(y ~ y') + G ( Ui , y|uj, y'), 



(All) 



L(y\y) = -j5(y - y') + G (y,y\y',y'). 



(A12) 
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APPENDIX B: TWO INTERACTING BOSONS IN ONE DIMENSION 



In this subsection an equation for the energy levels of a system with two interacting electrons in one dimension 
with periodic boundary conditions is derived. It turns out that this equation is the same as that obtained from the 
Bethe ansatz solution. Thus, the two methods lead to the same energy levels in this case. 

Consider a one dimensional box of length a, with periodic boundary conditions. The one-particle Hamiltonian is 
given by Eq. (||), the potential V(x) is zero. Using periodic boundary conditions the one-particle eigenfunctions of 
this Hamiltonian are given by 



and the eigenvalues are 



V fc (z) = -±= e 2mkx / a , (Bl) 



47T 2 fc 2 , . 

E k = (B2) 



where k — 0, ±1, ±2, • • •. From these wave functions one can construct the bosonic two-particle wave functions which 
may be written as 

f ^(*)tM»)+sMi/)tM*) for n>m 

i>n,m(x,y) = { ' (B3) 

{ ip n {x)ip m (y), for n = m. 

For two interacting particles the GF and the operator L{yi,y2) are given in Eqs. ( [l8|) and (po|), respectively. Using 
Eqs. (^) and (B3) the unperturbed t wo- particle GF and the operator L(yi, y 2 ) can be expressed in terms of the one 
particle eigenfunctions given in Eq. (Bl): 

n I n ,s _ 1 ^n{x)^ m {x')^* n {y)^*g{y') , (g) V>m (g0 V4 {vWn (j/) 

<^o(x,x \y,y ) - - p s p h - , {ai) 

z hi — h, n — hi m h — h, n — h/ m 

L(yuy 2 ) = -\ S( yi -y 2 )+X: (B5 ) 



where the strength of the interaction between the two particles is A, and E n is given in Eq. ( p32j ). 

The poles of the Green's function given in Eq. ( |Tg| ) are the energy eigenvalues of the system. The equation 
determining the poles of GF is 

det L( yi \y 2 ) = 0, (B6) 



where t he o perator L(yi\y2) is given by Eq. (B5). To evaluate the determinant it is convenient to use the basis given 
by Eq. (Bl). The matrix element Lki of the operator L(y±\y2) is then 

L ki = J o i>k(vi) L {yiiV2) Mv2) dy x d y2 = 8 kt i f-i + ^ ^2 E _ E ~ Ek — ^ ' ( B7 ) 

where we have used the integral 

1 

1>k( x ) V>nW ipm{x) dx = —= 5 n + m -k,o- (B8) 

V« 

Since the matrix Lki is diagonal, its determinant is 

+00 / +00 \ 

fe= — 00 \ ra= — 00 / 



+00 -, 

T= F T-2 " . (BIO) 

A ^ E - 14 (n 2 + (k - n) 2 ) V ' 



Thus, Eq. (p36|) is equivalent to 
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where k — 0, ±1, ±2, . . .. Using the following identities 



El 7T 7T 

z 2 - (2n)2 = 2l COt 2 Z ' 



n— — oo 



El 7T 7T 
— = ; -7T = tan— z. (Bll) 
z 2 - (2n + 1) 2 2z 2 v 7 



the right hand side of Eq. (BIO) can be written as 



y/2E- 


4ir 2 k 2 /a 2 


X 


V 2E - 


4:TT 2 k 2 /a 2 



2 - — = cot-y/2E- Air 2 k 2 /a 2 for k even, (B12) 



2 V t =-tan-v/2£-47r2fc2/ a 2 for odd. (B13) 

A 4 



This result is equivalent to the Bethe ansatz solution | ]19|J22| ] for a one-dimensional Bose gas with Dirac delta 
interactions. According to the Bethe ansatz solution, the energy eigenvalues E of the N interacting particles can be 
determined from the solution for kj of the following N equations: 



k t fe<j 



(-l) JV ~ 1 e"^ a = expU^-2arctan^-^ , (B14) 



A/2 

where j = 1, 2, . . . , N and 



iV 



E = J2tf- (B15) 

3=1 

For two particles, the two equations for k\ and ki can be rewritten as 

e- lPa = 1, (B16) 

e M ^°^ = exp f -4i arctan ^^^ j . (B17) 

where P = k\ + is the total momentum of the particles. From the first equation we have P = 2irn/a, where n is 
an integer, and the second equation may be written as a\j2E — P 2 + 2imi — —4 arctan v/2 f' / 2 ~i where n\ is another 
integer. Taking the tangent of this equation we find the same forms as given in Eqs. (B12)-(B13) for even and odd 
ri\, respectively. 
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